%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimation of a BVAR (draws for coefficient and VCV matrices)
% Minnesota and conjugate priors from Normal-Wishart
% based on code originally from: L.Gambetti and F. Canova, but changed and
% extended in a few ways
% this version: October 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inputs:
% y         = data
% nlagsvar  = Number of lags
% hndx      = All variables in the VAR are in first differences except for
% hours
% c         = 1: constant, 0: no constant
% ndraws    = Number of draws
% prior=1     Minnesota type with fixed hyperparameters
%             Normal prior for the VAR coefficients and fixed variance
% prior=2     Normal Inverted-Wishart, flat prior
% prior=3     Normal Inverted-Wishart, Minnesota structure as in Kadiyala and Karlsson
% T         = total number of lags after adjusting for laglength

% outputs:
% Sigma     draws for the posterior of the variance 
% Bet       draws from the posterior for beta in matrix form
% B and S   median values from the posterior distribution

function [Bet,Sigma,B,S]=bvarest(y,nlagsvar,hndx,c,ndraws,prior,T,phi,decay,do_level,trendbreak,break_choice)
% function [Bet,Sigma,B,S]=bvarest(y,nlagsvar,hndx,c,ndraws,prior,T,phi,decay,do_level,do_break,break_choice)

% new detrending
[t n]=size(y);
% if break_choice <= 3
%     dettrend = 0;
% elseif break_choice == 31
%     dettrend = 1;
% elseif break_choice == 32
%     dettrend = 2;
% elseif break_choice == 33
%     dettrend = 3;
% end
% [YY,xx,XX] = VARvariables_0711(y,c,trendbreak,nlagsvar,dettrend);

[YY,xx,XX] = VARvariables(y,c,nlagsvar);

Y=YY(1:T,:);
X=XX(1:T,:);
Bols=(X'*X)\(X'*Y);
res=Y-X*Bols;
Sols = (res'*res)/(T-1);

% decay:
hh = zeros(nlagsvar,1);
for lag=1:nlagsvar
    hh(lag) = lag^decay; % harmonic decay
%     hh(lag) = decay^(-lag+1); % geometric decay
end;

% Prior parameters
if prior == 1 
    % Prior for Beta
    % for XX data:
    B0 = zeros(c+nlagsvar*n,n);
    % all variables are estimated in first differences
    % hours may be specified in levels
    if do_level==1
        B0(c+hndx,hndx) = 1;
    end
    % vectorization
    B0 = reshape(B0,n*(c+nlagsvar*n),1);
    
    % Prior for Sigma
    % Sigma depends on tightness hyperparameters
    SB0=zeros(size(B0,1),size(B0,1));
    % scaling parameters as in Kadiyala and Karlsson:
    % residual standard error of univariate autoregressions
    s = zeros(n);
    for j=1:n
        xsep = xx(:,(j-1)*(c+nlagsvar)+1:c+j*nlagsvar);
        Bsep = (xsep'*xsep)\(xsep'*Y(:,j));
        ressep = Y(:,j)-xsep*Bsep;
        s(j) = std(ressep);
    end
    
    co=0;
    for i=1:n
        for l=1:nlagsvar;
            for j=1:n
                co=co+1;
                if (i==j)
                    SB0(co+c*i,co+c*i)=phi(1)/hh(l);
                else
                    SB0(co+c*i,co+c*i)=phi(1)*(phi(2)/hh(l))*(s(i)/s(j))^2;
                end
            end
        end
    end
    
    if c~=0
        for i=1:n
            dSB = diag(SB0);
            dSB((i-1)*(n*nlagsvar+c)+1:(i-1)*(n*nlagsvar+c)+c) = (s(1)^2)*phi(3);
            SB0 = diag(dSB);
        end;
    end
elseif prior == 2 
    % taken from Uhlig, JME 2005
    nu0=0; % degrees of freedom
    N0=zeros(size(X'*X,1),size(X'*X,2)); % variance factor
    S0=diag(ones(n,1)); 
    B0=zeros(n*nlagsvar+c,n); 
elseif prior==3
    % parameters according to Kadiyala and Karlsson
    nu0= n+2;
    
    % Omega in KK (Omega=inv(N))
    % needs to be such that together with S0 one obtains the Minnesota prior variance 
    % scaling parameter
    s = zeros(1,n);
    for j=1:n
        xsep = xx(:,(j-1)*(c+nlagsvar)+1:c+j*nlagsvar);
        Bsep = (xsep'*xsep)\(xsep'*Y(:,j));
        ressep = Y(:,j)-xsep*Bsep;
        s(j) = std(ressep);
    end
  
    N0=zeros(size(X'*X,1),size(X'*X,2));
    co=0;
    for l=1:nlagsvar;
        for j=1:n
            co=co+1;
            N0(co+c,co+c)=(phi(1)/hh(l))*(1/s(j))^2;
        end
    end
    
    if c~=0
        dN = diag(N0);
        dN(1:c) = phi(3); % Kadiyala and Karlsson value
        N0 = diag(dN);
    end;
    N0 = inv(N0);

    % prior for Sigma
    S0=diag(s')*(nu0-n-1);
    
    % prior for Beta
    B0 = zeros(c+nlagsvar*n,n);
    % all variables are estimated in first differences, except for hours
    if do_level==1
        B0(c+hndx,hndx) = 1;
    end
end

% Posterior parameters
if prior==1 
    SB=inv(inv(SB0)+kron(inv(Sols),X'*X));
    B=SB*(SB0\B0+kron(inv(Sols),X)'*vec(Y,2)); 
elseif prior==2 || prior==3
    % formulas for posterior parameters derived in Uhlig, 2005
    nu=nu0+t;
    N=N0+X'*X;
    B=N\(N0*B0+X'*X*Bols);
    S=nu0/nu*S0+t/nu*Sols+1/nu*(Bols-B0)'*N0*inv(N)*(X')*X*(Bols-B0);
    % in case of prior 2, B=Bols and S=Sols
end
        
% Generate draws from the posterior distribution
if prior==1
    vecB = zeros(size(B,1),size(B,2),ndraws);
    dummy = vec_to_mat(vecB(:,1),nlagsvar,n,c,1);
    Bet = zeros(size(dummy,1),size(dummy,2),ndraws);
    Sigma = zeros(size(Sols,1),size(Sols,2),ndraws);
    for d=1:ndraws
        % Draw VAR parameters from the posterior
        % coeffcients are normal
        vecB(:,d)=mvnrnd(B,SB,1)';
        Bet(:,:,d)=vec_to_mat(vecB(:,d),nlagsvar,n,c,1);
        Sigma(:,:,d)=Sols;
    end
    B=vec_to_mat(B,nlagsvar,n,c,1);
    S = Sols;
elseif prior==2 || prior==3
    dummy1 = vec(B,2);
    vecB = zeros(size(dummy1,1),size(dummy1,2),ndraws);
    dummy2 = vec_to_mat(vecB(:,1),nlagsvar,n,c,1);
    Bet = zeros(size(dummy2,1),size(dummy2,2),ndraws);
    Sigma = zeros(size(S,1),size(S,2),ndraws);
    for d=1:ndraws
        % draw for Sigma from Wishart
        Sigma(:,:,d)=wishrnd(inv(S)/nu,nu);
        Sigma(:,:,d)=inv(Sigma(:,:,d));      
        % Draw VAR parameters from normal using Sigma
        vecB(:,d)=mvnrnd(vec(B,2),kron(Sigma(:,:,d),inv(N)),1)';
        Bet(:,:,d)=vec_to_mat(vecB(:,d),nlagsvar,n,c,1);
    end
    B=B';
end